Introduction

This report describes the results of a preregistered study available at: https://osf.io/582wx.

Note also that this data has been cleaned beforehand. Three datasets were merged (joined) through an inner join—2 Qualtrics surveys and 1 Inquisit task—so as to keep only participants who at least participated at each step of the study. Missing data will be imputed later on. Duplicates were addressed with the rempsyc::best_duplicate function, which keeps the duplicate with the least amount of missing values, and in case of ties, takes the first occurrence.

Packages & Data

Packages

library(rempsyc)
library(dplyr)
library(interactions)
library(performance)
library(see)
library(report)
library(bestNormalize)
library(psych)
library(visdat)
library(missForest)
library(doParallel)

summary(report(sessionInfo()))

The analysis was done using the R Statistical language (v4.2.1; R Core Team, 2022) on Windows 10 x64, using the packages iterators (v1.0.14), doParallel (v1.0.17), interactions (v1.1.5), performance (v0.10.2), see (v0.7.4), report (v0.5.5.4), foreach (v1.5.2), bestNormalize (v1.8.3.9000), psych (v2.2.9), missForest (v1.5), rempsyc (v0.1.0.9), visdat (v0.5.3) and dplyr (v1.0.10).

Data

# Read data
data <- read.table("data/fulldataset.txt", sep = "\t", header = TRUE)

# Dummy-code group variable
data <- data %>% 
  mutate(condition_dum = ifelse(condition == "Mindfulness", 1, 0),
         condition = as.factor(condition))

cat(report_participants(data, threshold = 1))

379 participants (Mean age = 42.9, SD = 12.8, range: [21, 81]; Gender: 58.3% women, 40.1% men, 1.58% non-binary; Country: 100.00% United States of America; Race: 75.73% White, 10.55% Black or African American, 7.12% Asian, 4.22% Mixed, 2.37% other)

# Allocation ratio
report(data$condition)

x: 2 levels, namely Control (n = 191, 50.40%) and Mindfulness (n = 188, 49.60%)

Preparation

At this stage, we define a list of our relevant variables.

# Make list of DVs
col.list <- c("blastintensity", "blastduration", "blastintensity.duration",
              "KIMS", "BSCS", "BAQ")

Data cleaning

In this section, we are preparing the data for analysis: (a) taking care of preliminary exclusions, (b) checking for and exploring missing values, (d) imputing missing data with missForest, (e) computing scale means, and (f) extracting reliability indices for our scales.

Preliminary exclusions

First, we only want to keep those who agreed to keep their participation in the study after the debriefing.

data %>% 
  count(debrief.consent)
debrief.consent n
Yes, I accept to maintain my participation to this study. 379

Nobody to exclude based on consent.

Second, we know that we only want to keep participants who had at least an 80% success rate in the critical experimental manipulation task. Let’s see how many participants have less than an 80% success rate. Those with missing values for variable manipsuccessleft will also be excluded since they have not completed the critical experimental manipulation in this study.

data %>% 
    summarize(success.80 = sum(manipsuccessleft < .80, 
                               na.rm = TRUE),
              is.na = sum(is.na(manipsuccessleft)))
success.80 is.na
10 0

There’s 10 people with success smaller than 80%, let’s exclude them.

data <- data %>% 
    filter(manipsuccessleft >= .80)
cat(report_participants(data, threshold = 1))

369 participants (Mean age = 43.2, SD = 12.8, range: [21, 81]; Gender: 58.8% women, 39.6% men, 1.63% non-binary; Country: 100.00% United States of America; Race: 76.69% White, 9.49% Black or African American, 7.05% Asian, 4.34% Mixed, 2.44% other)

Let’s also exclude those who failed 2 or more attention checks (i.e., keep with those with a score of two or more).

data <- data %>% 
    mutate(att_check = rowSums(
      select(., att_check1, att_check2, att_check3)))

data %>% 
  count(att_check)
att_check n
1 5
2 17
3 347

There’s 5 more exclusions here.

data <- data %>% 
  filter(att_check >= 2)

cat(report_participants(data, threshold = 1))

364 participants (Mean age = 43.4, SD = 12.7, range: [21, 81]; Gender: 59.1% women, 39.3% men, 1.65% non-binary; Country: 100.00% United States of America; Race: 76.65% White, 9.62% Black or African American, 7.14% Asian, 4.40% Mixed, 2.20% other)

Explore missing data

Missing items

# Check for nice_na
nice_na(data, scales = c("BSCS", "BAQ", "KIMS"))
var items na cells na_percent na_max na_max_percent all_na
BSCS_1:BSCS_7 7 0 2548 0.00 0 0.00 0
BAQ_1:BAQ_12 12 0 4368 0.00 0 0.00 0
KIMS_1:KIMS_39 39 2 14196 0.01 2 5.13 0
Total 95 18 34580 0.05 2 2.11 0

No missing data for our scales of interest, yeah! Except 2 items for KIMS.

Patterns of missing data

Let’s check for patterns of missing data.

# Smaller subset of data for easier inspection
data %>%
  select(manualworkerId:att_check2_raw, 
         condition:condition_dum) %>%
  vis_miss

Little’s MCAR test

# Let's use Little's MCAR test to confirm
# We have to proceed by "scale" because the function can only
# support 30 variables max at a time
library(naniar)
data %>% 
  select(BSCS_1:BSCS_7) %>% 
  mcar_test
statistic df p.value missing.patterns
0 0 0 1 
# a p-value of 0 means the test failed because there's no missing values.

data %>% 
  select(BAQ_1:BAQ_12) %>% 
  mcar_test
statistic df p.value missing.patterns
0 0 0 1 
# a p-value of 0 means the test failed because there's no missing values.

data %>% 
  select(KIMS_1:KIMS_20) %>% 
  mcar_test
statistic df p.value missing.patterns
34.45728 19 0.0162216 2 
data %>% 
  select(KIMS_21:KIMS_39) %>% 
  mcar_test
statistic df p.value missing.patterns
18.04105 18 0.4529514 2

Impute missing data

…

Here, we impute missing data with the missForest package, as it is one of the best imputation methods.

Imputation

# Need character variables as factors
# "Error: Can not handle categorical predictors with more than 53 categories."
# So we have to temporarily remove IDs also...
new.data <- data %>% 
  select(-c(manualworkerId, embeddedworkerId, gender,
            att_check1, att_check2, att_check3)) %>% 
  mutate(across(where(is.character), as.factor))

# Parallel processing
registerDoParallel(cores = 4)

# Variables
set.seed(100)
data.imp <- missForest(new.data, verbose = TRUE, parallelize = "variables")
# Total time is 2 sec (4*0.5) - 4 cores

# Extract imputed dataset
new.data <- data.imp$ximp

There are some variables we don’t actually want to impute, like country. We want to keep those NAs in that case. Let’s add them back. We also want to add ID back.

# Add ID
new.data <- bind_cols(manualworkerId = data$manualworkerId, new.data)

# Add back the NAs in country, attention checks, etc.
data <- new.data %>% 
  mutate(country.ip = data$country.ip,
         gender = data$gender,
         att_check1 = data$att_check1, 
         att_check2 = data$att_check2,
         att_check3 = data$att_check3)

Details

Why impute the data? van Ginkel explains,

Regardless of the missingness mechanism, multiple imputation is always to be preferred over listwise deletion. Under MCAR it is preferred because it results in more statistical power, under MAR it is preferred because besides more power it will give unbiased results whereas listwise deletion may not, and under NMAR it is also the preferred method because it will give less biased results than listwise deletion.

van Ginkel, J. R., Linting, M., Rippe, R. C. A., & van der Voort, A. (2020). Rebutting existing misconceptions about multiple imputation as a method for handling missing data. Journal of Personality Assessment, 102(3), 297-308. https://doi.org/10.1080/00223891.2018.1530680

Why missForest? It outperforms other imputation methods, including the popular MICE (multiple imputation by chained equations). You also don’t end up with several datasets, which makes it easier for following analyses. Finally, it can be applied to mixed data types (missings in numeric & categorical variables).

Waljee, A. K., Mukherjee, A., Singal, A. G., Zhang, Y., Warren, J., Balis, U., … & Higgins, P. D. (2013). Comparison of imputation methods for missing laboratory data in medicine. BMJ open, 3(8), e002847. https://doi.org/10.1093/bioinformatics/btr597

Stekhoven, D. J., & Bühlmann, P. (2012). MissForest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1), 112-118. https://doi.org/10.1093/bioinformatics/btr597

Scale Means

…

Now that we have imputed the missing data, we are ready to calculate our scale means.

Self-control

# Reverse code items 2, 4, 6, 7
data <- data %>% 
  mutate(across(starts_with("BSCS"), .names = "{col}r"))
data <- data %>% 
  mutate(across(c(BSCS_2, BSCS_4, BSCS_6, BSCS_7), ~nice_reverse(.x, 5), .names = "{col}r"))

# Get mean BSCS
data <- data %>% 
  mutate(BSCS = rowMeans(select(., BSCS_1r:BSCS_7r)))

Trait aggression

# Reverse code item 7
data <- data %>% 
  mutate(across(starts_with("BAQ"), .names = "{col}r"))
data <- data %>% 
  mutate(across(BAQ_7, ~nice_reverse(.x, 7), .names = "{col}r"))

# Get sum of BAQ
data <- data %>% 
  mutate(BAQ = rowMeans(select(., BAQ_1r:BAQ_12r)))

Trait Mindfulness

# Reverse code items 3-4, 8, 11-12, 14, 16, 18, 20, 22, 23-24, 27-28, 31-32, 35-36
data <- data %>% 
  mutate(across(starts_with("KIMS"), .names = "{col}r"))
data <- data %>% 
  mutate(across(all_of(paste0("KIMS_", c(3:4, 8, 11:12, 14, 16, 18, 20,
                                         22:24, 27:28, 31:32, 35:36))), 
                ~nice_reverse(.x, 5), .names = "{col}r"))

# Get sum of KIMS
data <- data %>% 
  mutate(KIMS = rowMeans(select(., KIMS_1r:KIMS_39r)))

Intensity * Duration

# Create new variable blastintensity.duration
data <- data %>% 
  mutate(blastintensity.duration = blastintensity * blastduration)

Reliability

…

Now that we have reversed our items, we can get the alphas for our different scales.

Self-control

data %>% 
  select(BSCS_1r:BSCS_7r) %>% 
  omega()
## Loading required namespace: GPArotation

## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.84 
## G.6:                   0.84 
## Omega Hierarchical:    0.71 
## Omega H asymptotic:    0.8 
## Omega Total            0.88 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##            g   F1*   F2*   F3*   h2   u2   p2
## BSCS_1r 0.62  0.52             0.67 0.33 0.57
## BSCS_2r 0.68        0.49       0.70 0.30 0.65
## BSCS_3r 0.51  0.51             0.52 0.48 0.50
## BSCS_4r 0.65              0.39 0.59 0.41 0.73
## BSCS_5r 0.49  0.32        0.22 0.39 0.61 0.60
## BSCS_6r 0.66        0.45       0.64 0.36 0.68
## BSCS_7r 0.58              0.39 0.50 0.50 0.67
## 
## With Sums of squares  of:
##    g  F1*  F2*  F3* 
## 2.53 0.65 0.47 0.36 
## 
## general/max  3.9   max/min =   1.82
## mean percent general =  0.63    with sd =  0.08 and cv of  0.12 
## Explained Common Variance of the general factor =  0.63 
## 
## The degrees of freedom are 3  and the fit is  0.01 
## The number of observations was  364  with Chi Square =  3.83  with prob <  0.28
## The root mean square of the residuals is  0.01 
## The df corrected root mean square of the residuals is  0.02
## RMSEA index =  0.027  and the 10 % confidence intervals are  0 0.097
## BIC =  -13.87
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 14  and the fit is  0.43 
## The number of observations was  364  with Chi Square =  155.26  with prob <  0.000000000000000000000000064
## The root mean square of the residuals is  0.11 
## The df corrected root mean square of the residuals is  0.13 
## 
## RMSEA index =  0.166  and the 10 % confidence intervals are  0.144 0.191
## BIC =  72.7 
## 
## Measures of factor score adequacy             
##                                                  g   F1*   F2*   F3*
## Correlation of scores with factors            0.85  0.70  0.63  0.56
## Multiple R square of scores with factors      0.72  0.48  0.40  0.31
## Minimum correlation of factor score estimates 0.44 -0.03 -0.21 -0.37
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*  F2*  F3*
## Omega total for total scores and subscales    0.88 0.75 0.80 0.69
## Omega general for total scores and subscales  0.71 0.44 0.54 0.49
## Omega group for total scores and subscales    0.13 0.31 0.26 0.20

Trait aggression

data %>% 
  select(BAQ_1r:BAQ_12r) %>% 
  omega()

## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.82 
## G.6:                   0.86 
## Omega Hierarchical:    0.63 
## Omega H asymptotic:    0.71 
## Omega Total            0.88 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##             g   F1*   F2*   F3*   h2   u2   p2
## BAQ_1r   0.47        0.65       0.66 0.34 0.34
## BAQ_2r   0.43        0.57       0.51 0.49 0.35
## BAQ_3r   0.51        0.55       0.57 0.43 0.45
## BAQ_4r-             -0.43       0.21 0.79 0.00
## BAQ_5r   0.30        0.58       0.43 0.57 0.20
## BAQ_6r   0.49  0.33             0.39 0.61 0.62
## BAQ_7r   0.45  0.42             0.40 0.60 0.51
## BAQ_8r   0.63  0.53             0.68 0.32 0.58
## BAQ_9r   0.67  0.50             0.71 0.29 0.64
## BAQ_10r  0.51              0.50 0.51 0.49 0.51
## BAQ_11r  0.56              0.46 0.56 0.44 0.57
## BAQ_12r  0.47              0.44 0.45 0.55 0.50
## 
## With Sums of squares  of:
##    g  F1*  F2*  F3* 
## 2.84 0.87 1.69 0.68 
## 
## general/max  1.68   max/min =   2.51
## mean percent general =  0.44    with sd =  0.19 and cv of  0.43 
## Explained Common Variance of the general factor =  0.47 
## 
## The degrees of freedom are 33  and the fit is  0.46 
## The number of observations was  364  with Chi Square =  164.75  with prob <  0.0000000000000000002
## The root mean square of the residuals is  0.05 
## The df corrected root mean square of the residuals is  0.07
## RMSEA index =  0.105  and the 10 % confidence intervals are  0.089 0.121
## BIC =  -29.85
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 54  and the fit is  2.02 
## The number of observations was  364  with Chi Square =  721.93  with prob <  0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000014
## The root mean square of the residuals is  0.16 
## The df corrected root mean square of the residuals is  0.18 
## 
## RMSEA index =  0.184  and the 10 % confidence intervals are  0.173 0.197
## BIC =  403.48 
## 
## Measures of factor score adequacy             
##                                                  g   F1*  F2*   F3*
## Correlation of scores with factors            0.81  0.67 0.84  0.67
## Multiple R square of scores with factors      0.66  0.45 0.70  0.45
## Minimum correlation of factor score estimates 0.32 -0.10 0.40 -0.11
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*  F2*  F3*
## Omega total for total scores and subscales    0.88 0.81 0.69 0.74
## Omega general for total scores and subscales  0.63 0.50 0.30 0.41
## Omega group for total scores and subscales    0.18 0.31 0.39 0.33

Trait Mindfulness

data %>% 
  select(KIMS_1r:KIMS_39r) %>% 
  omega()

## Omega 
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip, 
##     digits = digits, title = title, sl = sl, labels = labels, 
##     plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option, 
##     covar = covar)
## Alpha:                 0.9 
## G.6:                   0.94 
## Omega Hierarchical:    0.35 
## Omega H asymptotic:    0.38 
## Omega Total            0.93 
## 
## Schmid Leiman Factor loadings greater than  0.2 
##              g   F1*   F2*   F3*   h2   u2   p2
## KIMS_1r   0.22  0.50             0.32 0.68 0.16
## KIMS_2r   0.33  0.65  0.29       0.61 0.39 0.18
## KIMS_3r   0.49        0.35  0.40 0.53 0.47 0.46
## KIMS_4r   0.27        0.68       0.54 0.46 0.14
## KIMS_5r   0.25  0.43             0.28 0.72 0.23
## KIMS_6r   0.36  0.63  0.23       0.59 0.41 0.22
## KIMS_7r   0.43       -0.22  0.47 0.47 0.53 0.39
## KIMS_8r-        0.34 -0.46       0.33 0.67 0.00
## KIMS_9r   0.20  0.45             0.27 0.73 0.14
## KIMS_10r  0.31  0.69  0.20       0.62 0.38 0.15
## KIMS_11r  0.33        0.25  0.29 0.26 0.74 0.43
## KIMS_12r  0.25        0.67       0.51 0.49 0.12
## KIMS_13r  0.26  0.52             0.35 0.65 0.19
## KIMS_14r  0.33  0.41  0.54       0.57 0.43 0.19
## KIMS_15r  0.38  0.28        0.28 0.30 0.70 0.49
## KIMS_16r  0.25        0.71       0.57 0.43 0.11
## KIMS_17r        0.58 -0.23       0.42 0.58 0.06
## KIMS_18r  0.32  0.43  0.53       0.57 0.43 0.19
## KIMS_19r  0.29       -0.32  0.34 0.33 0.67 0.26
## KIMS_20r       -0.27  0.53       0.37 0.63 0.04
## KIMS_21r  0.25  0.62             0.45 0.55 0.13
## KIMS_22r  0.24  0.27  0.59       0.47 0.53 0.12
## KIMS_23r  0.44        0.30  0.35 0.41 0.59 0.47
## KIMS_24r              0.53       0.33 0.67 0.07
## KIMS_25r        0.47             0.26 0.74 0.12
## KIMS_26r  0.28  0.49             0.34 0.66 0.23
## KIMS_27r  0.35              0.38 0.32 0.68 0.38
## KIMS_28r  0.21        0.72       0.56 0.44 0.07
## KIMS_29r  0.24  0.45             0.27 0.73 0.21
## KIMS_30r  0.22  0.48             0.30 0.70 0.16
## KIMS_31r  0.39        0.21  0.40 0.36 0.64 0.41
## KIMS_32r  0.22        0.68       0.52 0.48 0.10
## KIMS_33r  0.22  0.47             0.27 0.73 0.18
## KIMS_34r  0.26  0.61             0.44 0.56 0.15
## KIMS_35r  0.41        0.31  0.36 0.40 0.60 0.42
## KIMS_36r              0.65       0.45 0.55 0.07
## KIMS_37r  0.27  0.56 -0.20       0.44 0.56 0.17
## KIMS_38r  0.34  0.28 -0.26  0.31 0.36 0.64 0.32
## KIMS_39r  0.23  0.52             0.35 0.65 0.16
## 
## With Sums of squares  of:
##   g F1* F2* F3* 
## 3.2 5.9 5.6 1.4 
## 
## general/max  0.55   max/min =   4.2
## mean percent general =  0.21    with sd =  0.13 and cv of  0.63 
## Explained Common Variance of the general factor =  0.2 
## 
## The degrees of freedom are 627  and the fit is  5.33 
## The number of observations was  364  with Chi Square =  1850.88  with prob <  0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000048
## The root mean square of the residuals is  0.06 
## The df corrected root mean square of the residuals is  0.06
## RMSEA index =  0.073  and the 10 % confidence intervals are  0.069 0.077
## BIC =  -1846.64
## 
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 702  and the fit is  14.7 
## The number of observations was  364  with Chi Square =  5123.73  with prob <  0
## The root mean square of the residuals is  0.22 
## The df corrected root mean square of the residuals is  0.22 
## 
## RMSEA index =  0.132  and the 10 % confidence intervals are  0.128 0.135
## BIC =  983.92 
## 
## Measures of factor score adequacy             
##                                                   g  F1*  F2*   F3*
## Correlation of scores with factors             0.68 0.91 0.94  0.68
## Multiple R square of scores with factors       0.46 0.83 0.87  0.46
## Minimum correlation of factor score estimates -0.08 0.66 0.75 -0.08
## 
##  Total, General and Subset omega for each subset
##                                                  g  F1*  F2*  F3*
## Omega total for total scores and subscales    0.93 0.89 0.85 0.78
## Omega general for total scores and subscales  0.35 0.17 0.12 0.41
## Omega group for total scores and subscales    0.44 0.72 0.74 0.37

We are getting the “Some items ( KIMS_8r ) were negatively correlated with the total scale and probably should be reversed.” warning. Let’s check these items. In the raw data, item 17 is:

labels.part2$KIMS_8
  1. “I tend to evaluate whether my perceptions are right or wrong.”

Which corresponds to items 17 and 19 as described in the original paper (Baer, Smith, & Allen in 2004). They also make sense to not reverse-score, theoretically speaking; they seem to relate to mindfulness, not its absence. So it is kind of strange that psych::alpha is suggesting to reverse-score these items. Perhaps the scale is not meant to be used as a single factor since it contains four factor subscales.

t-tests

In this section, we will: (a) test assumptions of normality, (b) transform variables violating assumptions, (c) test assumptions of homoscedasticity, (d) identify and winsorize outliers, and (e) conduct the t-tests.

Normality

lapply(col.list, function(x) 
  nice_normality(data, 
                 variable = x, 
                 title = x,
                 group = "condition",
                 shapiro = TRUE,
                 histogram = TRUE))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 1 rows containing non-finite values (`stat_density()`).

## 
## [[5]]

## 
## [[6]]

Several variables are clearly skewed. Let’s apply transformations. But first, let’s deal with the working memory task, SOPT (Self-Ordered Pointing Task). It is clearly problematic.

Transformation

The function below transforms variables according to the best possible transformation (via the bestNormalize package), and also standardizes the variables.

predict_bestNormalize <- function(var) {
  x <- bestNormalize(var, standardize = FALSE, allow_orderNorm = FALSE)
  print(cur_column())
  print(x$chosen_transform)
  cat("\n")
  predict(x)
}

set.seed(100)
data <- data %>% 
  mutate(across(all_of(col.list), 
                predict_bestNormalize,
                .names = "{.col}.t"))
## [1] "blastintensity"
## I(x) Transformation with 364 nonmissing obs.
## 
## [1] "blastduration"
## I(x) Transformation with 364 nonmissing obs.
## 
## [1] "blastintensity.duration"
## Non-Standardized sqrt(x + a) Transformation with 364 nonmissing obs.:
##  Relevant statistics:
##  - a = 0 
##  - mean (before standardization) = 66.62117 
##  - sd (before standardization) = 32.39192 
## 
## [1] "KIMS"
## Non-Standardized asinh(x) Transformation with 363 nonmissing obs.:
##  Relevant statistics:
##  - mean (before standardization) = 1.93303 
##  - sd (before standardization) = 0.1366693 
## 
## [1] "BSCS"
## Non-Standardized Yeo-Johnson Transformation with 364 nonmissing obs.:
##  Estimated statistics:
##  - lambda = 1.59839 
##  - mean (before standardization) = 6.664788 
##  - sd (before standardization) = 1.998245 
## 
## [1] "BAQ"
## Non-Standardized sqrt(x + a) Transformation with 364 nonmissing obs.:
##  Relevant statistics:
##  - a = 0 
##  - mean (before standardization) = 1.712829 
##  - sd (before standardization) = 0.3038779
col.list <- paste0(col.list, ".t")

Note. The I(x) transformations above are actually not transformations, but a shorthand function for passing the data “as is”. Suggesting the package estimated the various attempted transformations did not improve normality in those cases, so no transformation is used. This only appears when standardize is set to FALSE. When set to TRUE, for those variables, it is actually center_scale(x), suggesting that the data are only CENTERED because they need no transformation (no need to be scaled), only to be centered.

Let’s check if normality was corrected.

# Group normality
lapply(col.list, function(x) 
  nice_normality(data, 
                 x, 
                 "condition",
                 shapiro = TRUE,
                 title = x,
                 histogram = TRUE))
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 1 rows containing non-finite values (`stat_density()`).

## 
## [[5]]

## 
## [[6]]

Looks rather reasonable now, though not perfect (fortunately t-tests are quite robust against violations of normality).

We can now resume with the next step: checking variance.

Homoscedasticity

# Plotting variance
plots(lapply(col.list, function(x) {
  nice_varplot(data, x, group = "condition")
  }),
  n_columns = 3)
## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

Variance looks good. No group has four times the variance of any other group. We can now resume with checking outliers.

Outliers

We check outliers visually with the plot_outliers function, which draws red lines at +/- 3 median absolute deviations.

plots(lapply(col.list, function(x) {
  plot_outliers(data, x, group = "condition", ytitle = x, binwidth = 0.15)
  }),
  n_columns = 2)
## Warning: Removed 1 rows containing missing values (`stat_bindot()`).

There are some outliers, but nothing unreasonable. Let’s still check with the 3 median absolute deviations (MAD) method.

data %>% 
  #as.data.frame %>% 
  filter(condition == "Control") %>% 
  find_mad(col.list, criteria = 3)
## There were no outlier based on 3 median absolute deviations.
data %>% 
  #as.data.frame %>% 
  filter(condition == "Mindfulness") %>% 
  find_mad(col.list, criteria = 3)
## 3 outlier(s) based on 3 median absolute deviations for variable(s): 
##  blastintensity.t, blastduration.t, blastintensity.duration.t, KIMS.t, BSCS.t, BAQ.t 
## 
## Outliers per variable: 
## 
## $KIMS.t
##   Row KIMS.t_mad
## 1 139  -4.610287
## 2 163  -3.030456
## 
## $BAQ.t
##   Row BAQ.t_mad
## 1 124  3.146018

There are 3 outliers after our transformations.

Multivariate outliers

For multivariate outliers, it is recommended to use the Minimum Covariance Determinant, a robust version of the Mahalanobis distance (MCD, Leys et al., 2019).

Leys, C., Delacre, M., Mora, Y. L., Lakens, D., & Ley, C. (2019). How to classify, detect, and manage univariate and multivariate outliers, with emphasis on pre-registration. International Review of Social Psychology, 32(1).

x <- check_outliers(na.omit(data[col.list]), method = "mcd")
x
## 74 outliers detected: cases 3, 9, 17, 24, 26, 27, 34, 37, 39, 43, 47,
##   54, 56, 68, 72, 74, 76, 78, 79, 80, 94, 98, 100, 102, 105, 112, 113,
##   125, 127, 128, 136, 142, 145, 151, 158, 161, 164, 168, 172, 177, 185,
##   189, 200, 204, 208, 210, 211, 216, 218, 229, 230, 234, 237, 246, 251,
##   272, 278, 283, 287, 289, 295, 296, 307, 316, 322, 323, 324, 328, 335,
##   336, 337, 340, 350, 360.
## - Based on the following method and threshold: mcd (22.46).
## - For variables: blastintensity.t, blastduration.t,
##   blastintensity.duration.t, KIMS.t, BSCS.t, BAQ.t.
plot(x)

There are 74 multivariate outliers according to the MCD method.

This time, instead of winsorizing outliers, we attempt to exclude them to see if it makes any difference. We only exclude multivariate outliers, not univariate ones.

data <- data[-which(x), ]

# Update col.list
data <- data %>% 
  mutate(across(all_of(col.list), 
                as.numeric,
                .names = "{.col}.w"))

col.list <- paste0(col.list, ".w")

Standardization

We can now standardize our variables.

data <- data %>%
  mutate(across(all_of(col.list),
                function(x) {
                  as.numeric(scale(x))
                  },
                .names = "{col}.s"))

# Update col.list
col.list <- paste0(col.list, ".s")

We are now ready to compare the group condition (Control vs. Mindfulness Priming) across our different variables with the t-tests.

t-tests

nice_t_test(data, 
            response = col.list, 
            group = "condition") %>% 
  nice_table(highlight = 0.10, width = .80)
## Using Welch t-test (base R's default; cf. https://doi.org/10.5334/irsp.82).
## For the Student t-test, use `var.equal = TRUE`. 
##  

Interpretation: There is no effect whatsoever on sound blast. However, there seems to be baseline group differences in terms of self-control, trait aggression, and trait mindfulness, despite the block randomization…

Violin plots

Intensity * Duration

nice_violin(data, 
            group = "condition", 
            response = "blastintensity.duration.t.w.s",
            comp1 = 1,
            comp2 = 2,
            obs = TRUE,
            has.d = TRUE,
            d.y = 1)

Blast Intensity

nice_violin(data, 
            group = "condition", 
            response = "blastintensity.t.w.s",
            comp1 = 1,
            comp2 = 2,
            obs = TRUE,
            has.d = TRUE,
            d.y = 1)

Blast Duration

nice_violin(data, 
            group = "condition", 
            response = "blastduration.t.w.s",
            comp1 = 1,
            comp2 = 2,
            obs = TRUE,
            has.d = TRUE,
            d.y = 1)

Means, SD

Let’s extract the means and standard deviations for journal reporting.

Intensity * Duration

data %>% 
    group_by(condition) %>% 
    summarize(M = mean(blastintensity.duration),
              SD = sd(blastintensity.duration),
              N = n()) %>% 
  nice_table(width = 0.40)

Blast Intensity

data %>% 
    group_by(condition) %>% 
    summarize(M = mean(blastintensity),
              SD = sd(blastintensity),
              N = n()) %>% 
  nice_table(width = 0.40)

Blast Duration

data %>% 
    group_by(condition) %>% 
    summarize(M = mean(blastduration),
              SD = sd(blastduration),
              N = n()) %>% 
  nice_table(width = 0.40)

Moderations (confirmatory)

Let’s see if our variables don’t interact together with our experimental condition. But first, let’s test the models assumptions.

Assumptions

Intensity * Duration

big.mod3 <- lm(blastintensity.duration.t.w.s ~ condition_dum*BSCS.t.w.s
                 # + condition_dum*KIMS.t.w.s + condition_dum*BAQ.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod3)
## Variable `Component` is not in your data frame :/

Blast Intensity

big.mod1 <- lm(blastintensity.t.w.s ~ condition_dum*BSCS.t.w.s
                 # condition_dum*KIMS.t.w.s + condition_dum*BAQ.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod1)
## Variable `Component` is not in your data frame :/

Blast Duration

big.mod2 <- lm(blastduration.t.w.s ~ condition_dum*BSCS.t.w.s
                 # + condition_dum*KIMS.t.w.s + condition_dum*BAQ.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod2)
## Variable `Component` is not in your data frame :/

All the models assumptions look pretty good overall actually, even with all these variables. The lines for linearity and homoscedasticity are a bit skewed but nothing too crazy. Let’s now look at the results.

Moderations

Intensity * Duration

big.mod3 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Intensity

big.mod1 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Duration

big.mod2 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Interpretation: The condition by trait self-control (brief self-control scale, BSCS) interaction does not come up.

Interaction plots

Let’s plot the main interaction(s).

Intensity * Duration

interact_plot(big.mod3, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Blast Intensity

interact_plot(big.mod1, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Blast Duration

interact_plot(big.mod2, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Interpretation: The interaction is pretty much the same for all models. Counterintuitively, for people with low self-control, the priming mindfulness condition relates to lower aggression relative to the control condition. In contrast, for people with high self-control, the priming mindfulness condition relates to higher aggression.

Simple slopes

Let’s look at the simple slopes now (only for the significant interaction).

Intensity * Duration

big.mod3 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Intensity

big.mod1 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Duration

big.mod2 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Interpretation: There seems to have no effect of priming mindfulness on blast intensity as a function of self-control.

Moderations (exploratory)

Let’s see if our variables don’t interact together with our experimental condition. But first, let’s test the models assumptions.

Assumptions

Intensity * Duration

big.mod3 <- lm(blastintensity.duration.t.w.s ~ condition_dum*KIMS.t.w.s*BSCS.t.w.s +
                 condition_dum*BAQ.t.w.s*BSCS.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod3)
## Variable `Component` is not in your data frame :/

Blast Intensity

big.mod1 <- lm(blastintensity.t.w.s ~ condition_dum*KIMS.t.w.s*BSCS.t.w.s + 
                 condition_dum*BAQ.t.w.s*BSCS.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod1)
## Variable `Component` is not in your data frame :/

Blast Duration

big.mod2 <- lm(blastduration.t.w.s ~ condition_dum*KIMS.t.w.s*BSCS.t.w.s + 
                 condition_dum*BAQ.t.w.s*BSCS.t.w.s
               , data = data, na.action="na.exclude")
check_model(big.mod2)
## Variable `Component` is not in your data frame :/

All the models assumptions look pretty good overall actually, even with all these variables. The lines for linearity and homoscedasticity are a bit skewed but nothing too crazy. Let’s now look at the results.

Moderations

Intensity * Duration

big.mod3 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Intensity

big.mod1 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Duration

big.mod2 %>% 
  nice_lm(b.label = "B") %>% 
  nice_table(highlight = TRUE)

Interpretation: The condition by trait self-control (brief self-control scale, BSCS) interaction does not come up.

Interaction plots

Let’s plot the main significant interaction(s).

Intensity * Duration

interact_plot(big.mod3, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Blast Intensity

interact_plot(big.mod1, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Blast Duration

interact_plot(big.mod2, pred = "condition_dum", modx = "BSCS.t.w.s", 
              modxvals = NULL, interval = TRUE, x.label = "Condition", 
              pred.labels = c("Control", "Mindfulness"),
              legend.main = "Trait Self-Control")

Interpretation: The interaction is pretty much the same for all models. Counterintuitively, for people with low self-control, the priming mindfulness condition relates to lower aggression relative to the control condition. In contrast, for people with high self-control, the priming mindfulness condition relates to higher aggression.

Simple slopes

Let’s look at the simple slopes now (only for the significant interaction).

Intensity * Duration

big.mod3 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Intensity

big.mod1 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Blast Duration

big.mod2 %>%
  nice_lm_slopes(predictor = "condition_dum",
                 moderator = "BSCS.t.w.s",
                 b.label = "B") %>% 
  nice_table(highlight = TRUE)

Interpretation: There seems to have no effect of priming mindfulness on blast intensity as a function of self-control.

Conclusions

Based on the results, it seems that the predicted interaction between self-control and the priming mindfulness manipulation does not come up.

Package References

report::cite_packages(sessionInfo())
LS0tDQp0aXRsZTogJyoqUHJpbWluZyBNaW5kZnVsbmVzcyBQcm9qZWN0IDIqKicNCnN1YnRpdGxlOiBDb21wYXJpc29uICYgYW5hbHlzaXMgcmVwb3J0DQphdXRob3I6ICJSw6ltaSBUaMOpcmlhdWx0Ig0KZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiDQpvdXRwdXQ6DQogIGh0bWxfZG9jdW1lbnQ6DQogICAgdGhlbWU6IGNlcnVsZWFuDQogICAgaGlnaGxpZ2h0OiBweWdtZW50cw0KICAgIHRvYzogeWVzDQogICAgdG9jX2RlcHRoOiAzDQogICAgdG9jX2Zsb2F0OiB5ZXMNCiAgICBudW1iZXJfc2VjdGlvbnM6IG5vDQogICAgZGZfcHJpbnQ6IGthYmxlDQogICAgY29kZV9mb2xkaW5nOiBzaG93ICMgb3I6IGhpZGUNCiAgICBjb2RlX2Rvd25sb2FkOiB5ZXMNCiAgICBhbmNob3Jfc2VjdGlvbnM6DQogICAgICBzdHlsZTogc3ltYm9sDQotLS0NCg0KYGBge3Igc2V0dXAsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9VFJVRSwgaW5jbHVkZT1GQUxTRX0NCmZhc3QgPC0gRkFMU0UgICMgTWFrZSB0aGlzIHRydWUgdG8gc2tpcCB0aGUgY2h1bmtzDQpgYGANCg0KYGBge3Iga2xpcHB5LCBlY2hvPUZBTFNFLCBpbmNsdWRlPVRSVUV9DQprbGlwcHk6OmtsaXBweShwb3NpdGlvbiA9IGMoJ3RvcCcsICdyaWdodCcpKQ0KYGBgDQoNCiMgSW50cm9kdWN0aW9uDQoNClRoaXMgcmVwb3J0IGRlc2NyaWJlcyB0aGUgcmVzdWx0cyBvZiBhIHByZXJlZ2lzdGVyZWQgc3R1ZHkgYXZhaWxhYmxlIGF0OiBodHRwczovL29zZi5pby81ODJ3eC4NCg0KLS0tDQpOb3RlIGFsc28gdGhhdCB0aGlzIGRhdGEgaGFzIGJlZW4gY2xlYW5lZCBiZWZvcmVoYW5kLiBUaHJlZSBkYXRhc2V0cyB3ZXJlIG1lcmdlZCAoam9pbmVkKSB0aHJvdWdoIGFuIGlubmVyIGpvaW7igJQyIFF1YWx0cmljcyBzdXJ2ZXlzIGFuZCAxIElucXVpc2l0IHRhc2vigJRzbyBhcyB0byBrZWVwIG9ubHkgcGFydGljaXBhbnRzIHdobyBhdCBsZWFzdCBwYXJ0aWNpcGF0ZWQgYXQgZWFjaCBzdGVwIG9mIHRoZSBzdHVkeS4gTWlzc2luZyBkYXRhIHdpbGwgYmUgaW1wdXRlZCBsYXRlciBvbi4gRHVwbGljYXRlcyB3ZXJlIGFkZHJlc3NlZCB3aXRoIHRoZSBgcmVtcHN5Yzo6YmVzdF9kdXBsaWNhdGVgIGZ1bmN0aW9uLCB3aGljaCBrZWVwcyB0aGUgZHVwbGljYXRlIHdpdGggdGhlIGxlYXN0IGFtb3VudCBvZiBtaXNzaW5nIHZhbHVlcywgYW5kIGluIGNhc2Ugb2YgdGllcywgdGFrZXMgdGhlIGZpcnN0IG9jY3VycmVuY2UuDQoNCg0KIyBQYWNrYWdlcyAmIERhdGENCg0KIyMgUGFja2FnZXMNCg0KYGBge3Igd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgcmVzdWx0cz0nYXNpcyd9DQpsaWJyYXJ5KHJlbXBzeWMpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShpbnRlcmFjdGlvbnMpDQpsaWJyYXJ5KHBlcmZvcm1hbmNlKQ0KbGlicmFyeShzZWUpDQpsaWJyYXJ5KHJlcG9ydCkNCmxpYnJhcnkoYmVzdE5vcm1hbGl6ZSkNCmxpYnJhcnkocHN5Y2gpDQpsaWJyYXJ5KHZpc2RhdCkNCmxpYnJhcnkobWlzc0ZvcmVzdCkNCmxpYnJhcnkoZG9QYXJhbGxlbCkNCg0Kc3VtbWFyeShyZXBvcnQoc2Vzc2lvbkluZm8oKSkpDQpgYGANCg0KIyMgRGF0YQ0KDQpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPVRSVUUsIHJlc3VsdHM9J2FzaXMnfQ0KIyBSZWFkIGRhdGENCmRhdGEgPC0gcmVhZC50YWJsZSgiZGF0YS9mdWxsZGF0YXNldC50eHQiLCBzZXAgPSAiXHQiLCBoZWFkZXIgPSBUUlVFKQ0KDQojIER1bW15LWNvZGUgZ3JvdXAgdmFyaWFibGUNCmRhdGEgPC0gZGF0YSAlPiUgDQogIG11dGF0ZShjb25kaXRpb25fZHVtID0gaWZlbHNlKGNvbmRpdGlvbiA9PSAiTWluZGZ1bG5lc3MiLCAxLCAwKSwNCiAgICAgICAgIGNvbmRpdGlvbiA9IGFzLmZhY3Rvcihjb25kaXRpb24pKQ0KDQpjYXQocmVwb3J0X3BhcnRpY2lwYW50cyhkYXRhLCB0aHJlc2hvbGQgPSAxKSkNCg0KIyBBbGxvY2F0aW9uIHJhdGlvDQpyZXBvcnQoZGF0YSRjb25kaXRpb24pDQoNCmBgYA0KDQojIyBQcmVwYXJhdGlvbg0KDQpBdCB0aGlzIHN0YWdlLCB3ZSBkZWZpbmUgYSBsaXN0IG9mIG91ciByZWxldmFudCB2YXJpYWJsZXMuDQoNCmBgYHtyIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9VFJVRSwgcmVzdWx0cz0nYXNpcyd9DQojIE1ha2UgbGlzdCBvZiBEVnMNCmNvbC5saXN0IDwtIGMoImJsYXN0aW50ZW5zaXR5IiwgImJsYXN0ZHVyYXRpb24iLCAiYmxhc3RpbnRlbnNpdHkuZHVyYXRpb24iLA0KICAgICAgICAgICAgICAiS0lNUyIsICJCU0NTIiwgIkJBUSIpDQoNCmBgYA0KDQojIERhdGEgY2xlYW5pbmcNCg0KYGBge3IgY2hpbGQ9aWYgKGZhc3QgPT0gRkFMU0UpICcwX2NsZWFuaW5nLlJtZCcsIGV2YWwgPSBUUlVFfQ0KYGBgDQoNCiMgdC10ZXN0cw0KDQpgYGB7ciBjaGlsZD1pZiAoZmFzdCA9PSBGQUxTRSkgJzFfdF90ZXN0cy5SbWQnLCBldmFsID0gVFJVRX0NCmBgYA0KDQojIE1vZGVyYXRpb25zIChjb25maXJtYXRvcnkpDQoNCmBgYHtyIGNoaWxkPWlmIChmYXN0ID09IEZBTFNFKSAnMl9tb2RlcmF0aW9ucy5SbWQnLCBldmFsID0gVFJVRX0NCmBgYA0KDQojIE1vZGVyYXRpb25zIChleHBsb3JhdG9yeSkNCg0KYGBge3IgY2hpbGQ9aWYgKGZhc3QgPT0gRkFMU0UpICczX21vZGVyYXRpb25zMi5SbWQnLCBldmFsID0gVFJVRX0NCmBgYA0KDQojIENvbmNsdXNpb25zDQoNCkJhc2VkIG9uIHRoZSByZXN1bHRzLCBpdCBzZWVtcyB0aGF0IHRoZSBwcmVkaWN0ZWQgaW50ZXJhY3Rpb24gYmV0d2VlbiBzZWxmLWNvbnRyb2wgYW5kIHRoZSBwcmltaW5nIG1pbmRmdWxuZXNzIG1hbmlwdWxhdGlvbiBkb2VzIG5vdCBjb21lIHVwLg0KDQojIFBhY2thZ2UgUmVmZXJlbmNlcw0KDQpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCByZXN1bHRzPSdhc2lzJ30NCnJlcG9ydDo6Y2l0ZV9wYWNrYWdlcyhzZXNzaW9uSW5mbygpKQ0KYGBgDQo=